function [x, U, times] = fisher()

%  [x,u,t,info] = burgers;

%FISHER Solve Fisher's equation using the trapezoidal rule in time.
%
%   [x,U,times] = FISHER() solves
%
%       u_t = ep*u_xx + u*(1-u)
%
%   on x0 < x < x1. The trapezoidal (Crank-Nicolson) rule is
%   used in time, and Newton's method is used at each time step.
%   Each Newton correction is obtained with a tridiagonal Thomas solve.
%
%   The columns of U contain the solution at the times in TIMES.
%   The program also writes the data to fisher.txt.

    ep     = 0.1;
    nx     = 2000;
    lambda = 1.0;
    x0     = -10.0;
    x1     = 20.0;
    nt     = 1000;
    tmax   = 10.0;

    dx = (x1 - x0)/(nx + 1);
    dt = tmax/nt;

    it1 = floor(nt/2);
    it2 = nt;

    % Interior grid points.
    x = x0 + (1:nx)'*dx;

    % Initial condition and its second derivative.
    ex  = exp(lambda*x);
    g   = 1./(1 + ex);
    gxx = -lambda^2 .* ex .* (1 - ex) .* g.^3;

    % u includes the two boundary values. In the Fortran code, the
    % interior initial values of u were not explicitly assigned; they
    % must be initialized to g for Newton's method to have a valid start.
    u = zeros(nx + 2,1);
    u(1)       = 1/(1 + exp(lambda*x0));
    u(2:nx+1)  = g;
    u(nx+2)    = 1/(1 + exp(lambda*x1));

    q = ep*gxx + g.*(1 - g + 2/dt);

    snapshots = zeros(nx,2);

    for it = 1:nt
        t = it*dt;

        arg0 = lambda*x0 + ep*t*lambda^2 - t;
        arg1 = lambda*x1 + ep*t*lambda^2 - t;
        u0 = 1/(1 + exp(arg0));
        u1 = 1/(1 + exp(arg1));

        u = solveBVP(nx,dx,dt,ep,u0,u1,u,q);

        if it == it1
            snapshots(:,1) = u(2:nx+1);
        end
        if it == it2
            snapshots(:,2) = u(2:nx+1);
        end

        uxx = (u(3:nx+2) - 2*u(2:nx+1) + u(1:nx))/dx^2;
        q = ep*uxx + u(2:nx+1).*(1 - u(2:nx+1) + 2/dt);
    end

    times = [0, it1*dt, it2*dt];
    U = [g, snapshots];

    fid = fopen('fisher.txt','w');
    if fid == -1
        error('fisher:FileOpenError', ...
              'Unable to open fisher.txt for writing.');
    end
    cleanupObject = onCleanup(@() fclose(fid)); 

    fprintf(fid,'     ep = %.4g     t1 = %10.4f     t2 = %10.4f\n', ...
            ep,times(2),times(3));
    fprintf(fid,'     x     t=0     t=t1     t=t2\n');
    fprintf(fid,'%14.4e     %14.4e     %14.4e     %14.4e\n', ...
            [x,U].');

    if nargout == 0
        plot(x,U,'LineWidth',1.2)
        xlabel('x')
        ylabel('u(x,t)')
        legend(compose('t = %g',times),'Location','best')
        grid on
    end
end

function y = solveBVP(nx,dx,dt,ep,y0,y1,y,q)
%SOLVEBVP Newton solve for one implicit trapezoidal time step.

    tolerance = 1.0e-5;
    maxIterations = 50;
    dx2 = dx^2;

    y(1) = y0;
    y(nx+2) = y1;

    lower = -ones(nx-1,1);
    upper = -ones(nx-1,1);

    for iteration = 1:maxIterations
        center = y(2:nx+1);
        left   = y(1:nx);
        right  = y(3:nx+2);

        f = -(center.*(1 - 2/dt - center) + q)/ep;
        residual = 2*center - right - left + dx2*f;

        fy = (-1 + 2/dt + 2*center)/ep;
        diagonal = 2 + dx2*fy;

        correction = tridiagonalSolve(lower,diagonal,upper,-residual);
        y(2:nx+1) = center + correction;

        if norm(correction,inf) < tolerance
            return
        end
    end

    warning('fisher:NewtonNotConverged', ...
            'Newton iteration did not converge in %d iterations.', ...
            maxIterations);
end

function x = tridiagonalSolve(lower,diagonal,upper,rhs)
%TRIDIAGONALSOLVE Solve a tridiagonal system by the Thomas algorithm.

    n = numel(diagonal);
    modifiedUpper = zeros(n-1,1);
    modifiedRhs   = zeros(n,1);

    pivot = diagonal(1);
    if pivot == 0
        error('fisher:ZeroPivot','Zero pivot in tridiagonal solve.');
    end

    if n > 1
        modifiedUpper(1) = upper(1)/pivot;
    end
    modifiedRhs(1) = rhs(1)/pivot;

    for k = 2:n
        pivot = diagonal(k) - lower(k-1)*modifiedUpper(k-1);
        if pivot == 0
            error('fisher:ZeroPivot','Zero pivot in tridiagonal solve.');
        end

        if k < n
            modifiedUpper(k) = upper(k)/pivot;
        end
        modifiedRhs(k) = ...
            (rhs(k) - lower(k-1)*modifiedRhs(k-1))/pivot;
    end

    x = zeros(n,1);
    x(n) = modifiedRhs(n);
    for k = n-1:-1:1
        x(k) = modifiedRhs(k) - modifiedUpper(k)*x(k+1);
    end
end
